Jean-Christophe Loiseau
Maître de Conférences
📧 : jean-christophe.loiseau@ensam.eu
📅 : December 19th 2019
# --> Julia packages for the plots.
using Plots
using LaTeXStrings
# --> Julia package to integrate forward in time ODEs.
using DifferentialEquations
# -->
using DynamicalSystems
The aim of this notebook, in association with the first lecture of the dynamical system class I teach at Arts et Métiers, is to illustrate some of the various concepts that will be taught during the course on a set of relatively simple, yet illustrative, examples. For that purpose, the following systems will be considered:
Double pendulum: canonical example from classical physics belonging to the class of Hamiltonian systems, i.e. systems for which a quantity (namely the Hamiltonian, here the sum of the potential and kinetic energies) is conserved.Lorenz system: canonical example from dynamical system theory belonging to the class of dissipative systems of particular historical importance.Logistic map: Another classical example from dynamical system theory used to illustrate the concept of bifurcations. Note that, contrary to all other examples considered herein, this one is formulated in a discrete-time framework.Provided you have installed julia and the required packages (i.e. Plots.jl, PyPlot.jl, DifferentialEquations.jl and DynamicalSystems.jl), this notebook should be self-contained.
All of the pieces of code needed to run the various examples are included.
If you encounter any problem when running this notebook, please do not hesitate to get back to me either by email (📧) or directly on GitHub and I'll try to fix it.
License: This work is distributed under the Creative Commons BY-NC-SA 4.0 license.
The double pendulum is a canonical example from classical mechanics which can be used to model various mechanical systems. It belongs to the class of so-called Hamiltonian systems for which a quantity (namely the Hamiltonian, here the sum of kinetic and potential energies) is conserved. Despite its apparent simplicity, this system can nonetheless exhibit extremely complex dynamics. Defining the angular positions $\theta_1$ and $\theta_2$ as shown in the left figure and the corresponding angular velocities $\dot{\theta}_1$ and $\dot{\theta}_2$, the Lagrangian $\mathcal{L}(\theta_1, \theta_2, \dot{\theta}_1, \dot{\theta}_2)$ of this system is the difference of the kinetic energy given by
$$ \mathcal{T} = \frac{1}{2} m_1 l_1^2 \dot{\theta}_1^2 + \frac{1}{2} m_2 \left( l_1^2 \dot{\theta}_1^2 + l_2^2 \dot{\theta}_2^2 + 2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \cos \left( \theta_1 - \theta_2 \right) \right) $$and potential energy
$$ \mathcal{V} = (m_1 + m_2) g l_1 \cos \theta_1 + m_2 g l_2 \cos \theta_2 $$such that
$$ \mathcal{L} = \mathcal{T} - \mathcal{V}. $$The equations of motion are then given by
$$ \frac{\mathrm{d}}{\mathrm{d}t} \left( \frac{\partial \mathcal{L}}{\partial \dot{\theta}_i} \right) - \frac{\partial \mathcal{L}}{\partial \theta_i} = 0. $$Similarly, the equations of motions could be derived from the Hamiltonian $\mathcal{H}(p_1, p_2, \theta_1, \theta_2)$
$$ \mathcal{H} = \frac{m_2 l_2^2 p_1^2 + (m_1 + m_2) p_2^2 - 2 m_2 l_1 l_2 p_1 p_2 \cos(\theta_1 - \theta_2)}{2 l_1^2 l_2^2 m_2 (m_1 + m_2 \sin^2(\theta_1 - \theta_2)} - (m_1 + m_2) g l_1 \cos \theta_1 - m_2 g l_2 \cos \theta_2 $$with
$$ p_i = \frac{\partial \mathcal{L}}{\partial \dot{\theta}_i} $$the generalized momenta. The equations of motion are then given by
$$ \begin{aligned} & \frac{\mathrm{d} \theta_i}{\mathrm{d}t} = \frac{\partial \mathcal{H}}{\partial p_i} \\ & \frac{\mathrm{d} p_i}{\mathrm{d}t} = - \frac{\partial \mathcal{H}}{\partial \theta_i}. \end{aligned} $$This Hamiltonian formulation of the problem will be used in julia below to simulate this system.
The next few cells set up the problem and simulate the evolution of the system for various initial conditions.
# --> Define the mass of each pendulum.
m₁ = 1.0
m₂ = 1.0
# --> Define the length of each pendulum.
l₁ = 1.0
l₂ = 1.0
# --> Intensity of the gravity.
g = 9.81
# --> Pack all the parameters in a tuple for the sake of simplicity.
params = (m₁, m₂, l₁, l₂, g)
function Hamiltonian_double_pendulum(p, q, params)
"""
This function computes the Hamiltonian for the double pendulum.
INPUTS
------
p :
q : Angular positions of the two pendulums.
params : Parameters defining the system.
OUTPUT
------
H : Scalar
Hamiltonian of the system.
"""
# --> Unpack parameters.
m₁, m₂, l₁, l₂, g = params
# --> Compute the Hamiltonian.
H₁ = m₂ * l₂^2 * p[1].^2 + (m₁ + m₂) * l₁^2 * p[2].^2 - 2m₂ * l₁ * l₂ * p[1] * p[2] * cos(q[1] - q[2])
H₂ = 2m₂ * l₁^2 * l₂^2 * (m₁ + m₂ * sin(q[1] - q[2])^2)
H₃ = -(m₁ + m₂) * g * l₁ * cos(q[1]) - m₂ * g * l₂ * cos(q[2])
H = H₁ / H₂ + H₃
return H
end
# --> Temporal window over which the system is simulated.
tspan = (0.0, 25.0)
# --> Instants of time at which the evolution is sampled.
Δt = 0.025;
t = collect(tspan[1]:Δt:tspan[2]);
# --> Zero initial condition for the angular velocities.
# NOTE: Same condition for all simulations.
p₀ = [0.0; 0.0]
# --> Small amplitudes simulation.
q₀ = [π/5 ; π/5]
prob₀ = HamiltonianProblem(
Hamiltonian_double_pendulum,
p₀, q₀,
tspan,
p=params
)
sol₀ = solve(prob₀, McAte5(), tol=1e-8, dt=1e-3);
# --> Intermediate amplitudes simulation.
q₁ = [2π/5 ; 2π/5]
prob₁ = HamiltonianProblem(
Hamiltonian_double_pendulum,
p₀, q₁,
tspan,
p=params
)
sol₁ = solve(prob₁, McAte5(), tol=1e-8, dt=1e-3);
# --> Intermediate amplitudes simulation (bis).
q₂ = [3π/5 ; 3π/5]
prob₂ = HamiltonianProblem(
Hamiltonian_double_pendulum,
p₀, q₂,
tspan,
p=params
)
sol₂ = solve(prob₂, McAte5(), tol=1e-8, dt=1e-3);
# --> Large amplitudes simulation
q₃ = [4π/5 ; 4π/5]
prob₃ = HamiltonianProblem(
Hamiltonian_double_pendulum,
p₀, q₃,
tspan,
p=params
)
sol₃ = solve(prob₃, McAte5(), tol=1e-8, dt=1e-3);
function polar_to_cartesian(t, sol, params)
"""
Utility function to go from polar to cartesian coordinates.
"""
# --> Unpack parameters.
m₁, m₂, l₁, l₂, g = params
# -->
x₁ = l₁ * sin.(sol(t)[3, :])
y₁ = -l₁ * cos.(sol(t)[3, :])
# -->
x₂ = x₁ + l₂ * sin.(sol(t)[4, :])
y₂ = y₁ - l₂ * cos.(sol(t)[4, :])
return [x₁ y₁ x₂ y₂]
end
# --> Transform each simulation into the cartesian reference frame.
x₀ = polar_to_cartesian(t, sol₀, params)
x₁ = polar_to_cartesian(t, sol₁, params)
x₂ = polar_to_cartesian(t, sol₂, params)
x₃ = polar_to_cartesian(t, sol₃, params);
function plot_pendulum(i, x)
"""
Utility function to plot the evolution of the double pendulum.
"""
# --> Plot the position of the first pendulum.
p = plot(
[0, x[i, 1]], [0, x[i, 2]],
size=(512, 512),
xlim=(-2.2, 2.2),
ylim=(-2.2, 2.2),
markersize=15,
markershape=:circle,
label="",
framestyle=:none,
)
# --> Plot the position of the second pendulum.
p = plot!(
[x[i, 1], x[i, 3]], [x[i, 2], x[i, 4]],
markersize=15,
markershape=:circle,
label="",
aspect_ratio=:equal,
)
# --> Plot the streak line effect.
if i > 9*2
p = plot!([x[i-3*2:i, 3]], [x[i-3*2:i, 4]],alpha = 0.15,linewidth = 2, color = :red,label ="");
p = plot!([x[i-5*2:i-3*2, 3]], [x[i-5*2:i-3*2, 4]],alpha = 0.08,linewidth = 2, color = :red,label ="");
p = plot!([x[i-7*2:i-5*2, 3]], [x[i-7*2:i-5*2, 4]],alpha = 0.04,linewidth = 2, color = :red, label ="");
p = plot!([x[i-9*2:i-7*2, 3]], [x[i-9*2:i-7*2, 4]],alpha = 0.01,linewidth = 2, color = :red, label="");
end
return p
end
gr()
# --> Loop through the frames.
anim = @animate for i = 1:length(t)
# --> Plot each simulation.
p₀ = plot_pendulum(i, x₀);
p₁ = plot_pendulum(i, x₁);
p₂ = plot_pendulum(i, x₂);
p₃ = plot_pendulum(i, x₃);
p = plot(
p₀, p₁, p₂, p₃,
layout=(2, 2),
size=(1024, 1024)
);
end
# --> Display the gif animation.
gif(anim, fps=30)